In this section, we show the implementation of a PDE-based model for pattern formation in Dicty, which was originally developed by Kessler & Levine, (1993) PRE,48(6):4801. In this model, cells are randomly placed on a 2D square lattice with density \(\rho\). Cyclic AMP concentration \(c\) (referred to as “con” in the code) obeys a reaction diffusion equation
\[\frac{\partial c}{\partial t} = a^2 (\frac{\partial^2 c}{\partial x^2}+\frac{\partial^2 c}{\partial y^2}) - kc + (\text{sources}) {\tag 1}\] where \(a\) is the lattice size, \(k\) is the degradation rate of cAMP.
Each cell has the capacity to alter its state as follows: it remains in an inactive state (state 1) until the concentration of cAMP surpasses a predefined threshold, denoted as \(c_T\) (referred to as “con_t” in the code). Once the cAMP concentration exceeds \(c_T\), the cell transitions to an excited state (state 2) and releases an amount \(\Delta c\) of cAMP (referred to as “dc” in the code) over a duration of \(t_e\) time units. Subsequently, after a period of \(t_e\), the cell enters a refractory state (state 3) for a duration of \(t_r\) time units, before eventually reverting to the inactive state (state 1). The term “sources” in Equation (1) describes the secretion of cAMP in state 2.
For the initial condition, we consider a linear gradient of cAMP along the y-axis, with an average concentration of 1. Cells occupy state 1 from one end of the x direction, whereas at the opposite end, a fraction \(p\) of cells are in state 3. The fraction of cells in state 3 increases linearly along the x-axis. We enforce a no-flux boundary condition in modeling diffusion of cAMP.
# model parameters
n = 101 # number of grid points per dimension;
rho = 0.2 # cell density
con_t = 1 # threshold concentration
dc = 300 # cAMP production during excited state
k = 0.5 # degradation rate
t_e = 2 # excitation time
t_r = 20 # refractory time
p = 1 # maximum fraction of refractory cells
t_total = 100 # total simulation time
dt = 0.01 # time step size
nframe = 100 # animation frames
set.seed(1) # random number seed
# Initialize cells
init_cells <- function(n, ncell, p, t_r){
# Input:
# n: number of grid points in each dimension
# ncell: number of cells
# p: maximum fraction of refractory cells
# t_r: refractory time
# Output:
# states: cell states, a vector of size ncell
# 1: inactive, con < con_th; 2: excited, con > con_th; 3: refractory
# --- cell excited for tau, then refractory for tr, before transit to 0
# ind : indices for x and y grid points for each cell, a matrix of (ncell, 2)
# con: concentration of each grid point, a matrix of (n, n)
# clocks: time count down for cell state transitions, a vector of size ncell
# Initialize cells
ind_cell = sort(sample(1:(n*n), ncell, replace = FALSE))
ind_cell_xy = matrix(0, nrow = ncell, ncol = 2)
j = 0
for(i in ind_cell){
j = j + 1
ind_x = (i-1)%%n+1
ind_y = as.integer((i-1)/n) + 1
ind_cell_xy[j,1] = ind_x
ind_cell_xy[j,2] = ind_y
}
# initial concentration
con_grad = seq(0,2,2/(n-1))
con = matrix(rep(con_grad,n), ncol = n, byrow = T)
# Assign initial cell states, random assignment
states = rep(1, ncell)
clocks = integer(ncell)
for(i in seq_len(ncell)){
ind_x = ind_cell_xy[i,1]
ind_y = ind_cell_xy[i,2]
if(runif(1) < ind_x/n*p){
states[i] = 3
clocks[i] = t_r # state 3 with t_r waiting time
}
}
return(list(states = states, ind = ind_cell_xy, con = con, clocks = clocks))
}
ncell = as.integer(n*n*rho)
dat = init_cells(n,ncell,p,t_r)
plot(dat$ind[,1], dat$ind[,2], col = dat$states, xlab = "x", ylab = "y", cex = 0.5, pch = 15)
# update cell states
update_states <- function(n, ncell, dat, con_t, t_e, t_r, dt){
for (i in seq_len(ncell)){
s = dat$states[i]
ind_x = dat$ind[i,1]
ind_y = dat$ind[i,2]
clock_now = dat$clocks[i]
if(s == 1){
if(dat$con[ind_x, ind_y] >= con_t){ # enter excited state
dat$states[i] = 2
dat$clocks[i] = t_e
}
}else if(s == 2){
if(clock_now >= dt){ # in excited state
dat$clocks[i] = clock_now - dt
}else{ # exit excited state, enter refractory state
dat$states[i] = 3
dat$clocks[i] = t_r
}
}else if(s == 3){
if(clock_now >= dt){ # in refractory state
dat$clocks[i] = clock_now - dt
}else{ # exit refractory state, enter inactive state
dat$states[i] = 1
dat$clocks[i] = 0
}
}
}
return(dat)
}
# main modeling function
dicty_modeling <- function(n, ncell, dat, con_t, dc, k, t_e, t_r, t_total, dt, nframe){
# n: number of grid points in each dimension
# ncell: number of cells
# dat: states (ncell), ind (ncell, 2), con (n, n), clocks(ncell)
# con_t: cAMP threshold concentration
# dc: cAMP production
# k: degradation rate
# t_e: excitation time
# t_r: refractory time
# t_total: total simulation time
# dt: time step size
# nframe: number of frames to save
nt_all = as.integer(t_total/dt)
rate_c = dc/t_e
t_ind_save = as.integer(nt_all/(nframe-1))
con = dat$con
l = 0
pos_cells = dat$ind
states_all = matrix(0, nframe, ncell)
states_all[l,] = dat$states
for (t_ind in seq_len(nt_all)){
dat = update_states(n, ncell, dat, con_t, t_e, t_r, dt)
c_x_plus_one = rbind(con[-1,],con[n,]) # con(X+dX,Y), no flux BC
c_x_minus_one = rbind(con[1,],con[-n,]) # con(X-dX,Y), no flux BC
c_y_plus_one = cbind(con[,-1],con[,n]) # con(X,Y+dY), no flux BC
c_y_minus_one = cbind(con[,1],con[,-n]) # con(X,Y-dY), no flux BC
dcdt = (c_x_plus_one + c_x_minus_one + c_y_plus_one + c_y_minus_one - 4*con) - k * con
for (i in seq_len(ncell)){
if(dat$states[i] == 2){
ind_x = dat$ind[i,1]
ind_y = dat$ind[i,2]
dcdt[ind_x,ind_y] = dcdt[ind_x, ind_y] + rate_c
}
}
con = con + dcdt * dt
dat$con = con
if(t_ind%%t_ind_save == 0){
l = l + 1
states_all[l,] = dat$states
}
}
return(states_all)
}
state_all_saved = dicty_modeling(n, ncell, dat, con_t, dc, k, t_e, t_r, t_total, dt, nframe)
frame = 10
plot(dat$ind[,1], dat$ind[,2], col = state_all_saved[frame,], xlab = "x", ylab = "y", cex = 0.5,pch = 15)
The following scripts iteratively generate plots in png format and convert them to an animation gif using ImageMagick.
library(png)
library(magick)
library(magrittr)
for(i in 1:nframe){
png(filename=sprintf("output_dicty_%05d.png",i), width = 960, height = 960)
plot(dat$ind[,1], dat$ind[,2], col = state_all_saved[i,], xlab = "x", ylab = "y", cex = 1.6, pch = 15)
dev.off()
}
list.files(path='.', pattern = 'output_dicty', full.names = TRUE) %>%
image_read() %>%
image_join() %>%
image_animate(fps=10) %>%
image_write("extra/data/07D/dicty.gif")
list.files(path='.', pattern = 'output_dicty', full.names = TRUE) %>%
file.remove()
Here’s the output animation showing spiral cell state changes